DALS018-统计模型1-二项分布与泊松分布，MLE

发表于 | 分类于
| 字数统计: 3,749 | 阅读时长 ≈ 15

前言

这一部分是《Data Analysis for the life sciences》的第7章统计模型的第1小节，这一部分的主要内容涉及高维数据统计的一些原理，例如二项分布，泊松分布，最大似然估计等。

相应的R markdown文档可以参考作者的Github。

当我们在文献中看到p值的时候，它意味着使用某种概率分布来对零假设进行量化。很多时候我们很容易地就能找到使用哪种分布。例如，在女士品茶这个案例中，我们就能使用简单的概率计睡来确定零分布。文献中的大多数p值都是基于线性模型的样本均值或最小二乘估计来进行计算的，并利用CLT近似地救出它们的零分布。

CLT有着理论支持，它能保证这种近似是正确的。但是，我们无法一直使用这种近似法，例如当我们的样本太少时。以前面部分中，我们提到了，当总体数据近似服从正态分布时，样本的平均近似服从t分布。但是这种假设没有理论上的支持。在这一部分里，我们就会对这种情况进行建模。在研究身高时，我们从经验中就知道，这是一个非常好的模型。

但是，这并不意味着我们收集的每个数据集都服从正态分布。我们常见的一些例子，例如掷硬币，中彩票的人数，美国的收入。正态分布不是建模过程中的唯一参数分布。在这一部分里，我们会描述一些最广泛使用的参数分布以及它们在生命科学研究中的作用。我们还会介绍贝叶斯统计相关的知识，以及给出使用分层模型的使用案例。

二项分布

在$N$次实验中，成功$S=k$次的概率，公式如下所示：

其中$p$是成功的概率，二项分布的一个最有名案例就是掷硬币，当我们掷$N$次硬币，有$S$次是正面朝上的概率，在这个案例中，$p=0.5$。

$S/N$是独立随机变量的均值，因此CLT告诉我们，当$N$足够大时，$S$近似地服从正态分布。这种分布在生命科学研究中应用很广，最近，此分布在NGS检测variant callers和表型(genotyper)应用很广，它的一个特殊分布则是泊松分布(poisson distribution)（二顶分布的极限是泊松分布）。

泊松分布

由于买彩票的结果只有两个，中彩票与不中彩票，因此，那些赢得彩票的人数从理论上服从二项分布（我们是假设每个人只能买一张彩票）的。试验次数$N$是买彩票的人数，这个数字通常非常大。但是，赢得彩票的的人数通常是0到3之间，这就表明，赢得彩票的人数是不服从正态分布的。那么，为什么CLT在这种情况下不成立呢？从数学上可以对此进行解释，但是直觉又告诉我们，如果成功的总和是如此接近，并且大于1，这种分布就不可能是正态分布，这里我们进行一个快速的模拟：

1
2
3
4
5
p=10^-7 ##1 in 10,000,0000 chances of winning
N=5*10^6 ##5,000,000 tickets bought
winners=rbinom(1000,N,p) ##1000 is the number of different lotto draws
tab=table(winners)
plot(tab)

1
2
3
4
> prop.table(tab)
winners
    0     1     2     3 
0.617 0.294 0.080 0.009

达种情况下，$N$会非常大，$p$会非常小，因此我们可以计算出$N \times p$（此时我们之为$\lambda$）的值，例如位于0到10之间，而$S$此时就服从泊松分布，其简单的公式如下所示：

RNA-seq分析中常常用到泊松分布。因为我们对数千个分子进行采样，并且多数基因仅仅代表了所有基因中一小部分，因此在RNA-seq分析中使用泊松分布比较合适。

我们知道了泊松分布后，它对我们有什么作用？其中一个方面就是这种分布可以为我们提供实际分析过程中关于统计属性的总结信息。例如我们只从患者组和对照组中各取一个样本进行RNA-seq实验并研究基因的倍数变化。此时，在零假设成立的前提下（也就是这两个样本没有差异），利用泊松分布，我们可以计算出，统计学上的统计变异取决于基因的总丰度。我们可以从数学上来展示一下这个过程，下面是模拟这个过程：

1
2
3
4
5
6
7
N=10000##number of genes
lambdas=2^seq(1,16,len=N) ##these are the true abundances of genes
y=rpois(N,lambdas)##note that the null hypothesis is true for all genes
x=rpois(N,lambdas)
ind=which(y>0 & x>0)##make sure no 0s due to ratio and log
library(rafalib)
splot(log2(lambdas),log2(y/x),subset=ind)

对于低较的lambda的值，存在着更多的变异，如果我们计算出了大于2倍或更高倍数的基因，那些对于低丰度的基因来说，假阳性率也会变大。

NGS与Poisson分布

在这一部分里，我们还是要使用公共数据，如下所示：

1
2
3
4
# BiocManager::install("parathyroidSE")
library(parathyroidSE) ##available from Bioconductor
data(parathyroidGenesSE)
se <- parathyroidGenesSE

上述的这个数据包含在SummarizedExperiment对象中，在此不再赘述，如下所示：

1
2
3
4
> class(parathyroidGenesSE)
[1] "RangedSummarizedExperiment"
attr(,"package")
[1] "SummarizedExperiment"

我们只需要知道它包含一个数据矩阵，其中每一行都是基因组特征（可以理解为基因），每列都是一个样本。我们可以使用assay()函数提取这个矩阵。对于这个数据集，数据矩阵中每个单元格的数值对应了一个特定样本的一个特定基因的reads数。因此，我们可以绘制出前面类似的图形，从而表明利用实验数据构建的模型的预测行为，如下所示：

1
2
3
4
x <- assay(se)[,23]
y <- assay(se)[,24]
ind=which(y>0 & x>0)##make sure no 0s due to ratio and log
splot((log2(x)+log2(y))/2,log(x/y),subset=ind)

如果我们计算出4个个体的标准差，它要比利用泊松模型预测的要高得多。假设大多数基因在不同的个体之间的表达不同，那么，如果泊松模型适用于此情况，那么此图中应该存在着线性关系：

1
2
3
4
5
6
library(rafalib)
library(matrixStats)
vars=rowVars(assay(se)[,c(2,8,16,21)]) ##we now these four are 4
means=rowMeans(assay(se)[,c(2,8,16,21)]) ##different individulsa
splot(means,vars,log="xy",subset=which(means>0&vars>0)) ##plot a subset of data
abline(0,1,col=2,lwd=2)

上面图中的变异包括生物变异，但是泊松分布则无法对此进行建模。此时，我们需要另外一种分布，即负二项分布（negative binomial distribution)，这种分布整合了服从泊松分布的采样变异(sampling variability)和生物变异(biological variability)，它更适合对这种情况进行建模。负二项分布有两个参数，在处理计数数据(count data)方面有更大的灵活性。有关RNA-seq与负二项分布的内容可以参考这篇文献《Differential expression analysis for sequence count data》。而泊松分布只是负二项分布的一个特例。

最大似然估计

为了说明最大似然估计(Maximum Likelihood Estimation, MLE)的概念，我们使用一个相对简单的数据集来进行演示，这个数据集包含HMCV（人巨细胞病毒，human cytomegalovirus）基因组的回文位置(palindrome location)信息。我们从HMCV基因组每4kb个间隔中读取回文序列的数目，如下所示：

1
2
3
4
5
6
7
8
datadir="http://www.biostat.jhsph.edu/bstcourse/bio751/data"
x=read.csv(file.path(datadir,"hcmv.csv"))[,2]
breaks=seq(0,4000*round(max(x)/4000),4000)
tmp=cut(x,breaks)
counts=table(tmp)
library(rafalib)
mypar(1,1)
hist(counts)

上图中的计数数据似乎服从泊松分布，但是，$\lambda$是多少呢？常见的估计$\lambda$的方法是最大似然估计(maximum likelihoood estimation)。最了找到最大的MLE，我们要注意这些数据是独立的，我们所观察到这些数据的概率可以用以下公式来表示：

当下面这个公式取最大值时候，MLE就等于$\lambda$：

在实际计算过程中我们通常使用log转换后似然值(log-likehood)。下面我们使用一些代码来计算一下任意的$\lambda$下的log-likehood，并且我们使用函数optimized()来计睡一下能够值得这个函数值最大的时候的$\lambda$，最后我们绘制出log-likehood的曲线，如下所示：

1
2
3
4
5
6
l<-function(lambda) sum(dpois(counts,lambda,log=TRUE))
lambdas<-seq(3,7,len=100)
ls <- exp(sapply(lambdas,l))
plot(lambdas,ls,type="l")
mle=optimize(l,c(0,10),maximum=TRUE)
abline(v=mle$maximum)

如果进行一些微积分计算，就会发现计算MLE其实非常简单，就是计算其均值，如下所示：

1
print( c(mle$maximum, mean(counts) ) )

结果如下所示：

1
2
> print( c(mle$maximum, mean(counts) ) )
[1] 5.157894 5.157895

我们会注意到，观察到的计数数据与通过泊松分布预测的计算数据非常吻合：

1
2
3
theoretical<-qpois((seq(0,99)+0.5)/100,mean(counts))
qqplot(theoretical,counts)
abline(0,1)

因此，我们可以使用泊松分布来计算回文数据，其中$\lambda = 5.16$。

正连续值的分布

不同的基因在生物学重复中的变化不同。在后面我们会介绍次聚类模型(hierarchical model)，这是在分析基因组学数学中最有影响力的统计学方法之一。该方法以能够极大地改善了分析差异基因的初级方法。它是通过模拟基因变异的分布来进行建模的。在这里我们会介绍这种方法的参数模型。

我们想模拟一下基因特异性标准误(gene-specific standard erros)。这些分布是正态分布吗？这里我们需要记住，我们现在模拟的是总体的标准误，因此CLT在此种情况下并不适用，即使我们我们拥有数几千个基因。

下面就是案例，在这个案例中，我们会使用一个实验数据，这个数据中包含了小鼠基因表达的技术重复和生物学重复，我们会计算出技术重复和生物学重复的基因特异性样本标准误：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
library(Biobase) ##available from Bioconductor
# install_github("genomicsclass/maPooling")
library(maPooling) ##available from course github repo
data(maPooling)
pd=pData(maPooling)
##determin which samples are bio reps and which are tech reps
strain=factor(as.numeric(grepl("b",rownames(pd))))
pooled=which(rowSums(pd)==12 & strain==1)
techreps=exprs(maPooling[,pooled])
individuals=which(rowSums(pd)==1 & strain==1)
##remove replicates
individuals=individuals[-grep("tr",names(individuals))]
bioreps=exprs(maPooling)[,individuals]
###now compute the gene specific standard deviations
library(matrixStats)
techsds=rowSds(techreps)
biosds=rowSds(bioreps)

下面提研究样本的标准差(sample standard deviation)：

1
2
3
4
5
6
###now plot
library(rafalib)
mypar()
shist(biosds,unit=0.1,col=1,xlim=c(0,1.5))
shist(techsds,unit=0.1,col=2,add=TRUE)
legend("topright",c("Biological","Technical"), col=c(1,2),lty=c(1,1))

从上面的图形中我们可以发现，生物学变异(biological variability)远远大于技术变异(technical variability)。这就为我们提供了一个证据，即基因确实具有基因特异性生物学变异。如果我们想要模拟这种变异，我们首先要注意到，这些变异不服从正态分布，因为上面的曲线明显出现了拖尾现象。此外，由于SD为正数，因此上面的分布也存在着一定的限制。现在我们使用qq图来看一下：

1
2
qqnorm(biosds)
qqline(biosds)

一些参数分布可以处理这些情况（例如严格正向右拖尾分布，strictly positive and heavy right tails）。其中比较典型的分布就是gamma分布与F分布，gamma分布的密度函数为：

这个分布有2个参数，分别是$\alpha$和$\beta$，它们间接控制了曲线的位置(location)与缩放(scale)。它们也控制着曲线的形状，有关这些分布的知识可以参考《Mathematical Statistics and Data Analysis》这本书。

gamma分布的两个特例是卡方分布与指数分布，我们前面了解到了可以使用卡方分布来计算二联表。对于卡方分布来说，我们的参数是$\alpha=\nu/2$，$\beta=2$，其中$\nu$表示自由度。对于指数分布来说，$\alpha=1$，$\beta=\lambda$。

F分布在ANOVA中使用，它也是一个正向分布，并且存在着严格的右拖尾曲线，它的2个参数影响其曲线形状，，F分布的密度函数为：

这个函数中含有$B$函数（即$\beta$函数），还有2个自由度，即$d_{1}$和$d_{2}$。

模拟变异

在后面部分中，我们会介绍用于改善变异估计的层级模型(hierarchical model)。在这些案例中，从数学上很方便地能模拟出$\sigma^2$的分布，关于层级模型的内容可以参考《Linear models and empirical bayes methods for assessing differential expression in microarray experiments.》这篇文献，此文献指出，对基因的标准差进行采样服从校正后的F统计：

其中$d$表示了计算中$s^2$的自由度。例如，当我们比较两组样本时，每组3个样本，那么自由度就是4（文中后面提到了This leaves two free parameters to adjust to the data，直译就是剩下的2个自由参数用于调整数据，不太理解）。这里的$d$用于控制位置(location)，$s_{0}$用于控制缩放(scale)。下面的代码是F分布的一些案例，我们绘制出了相应的图形，并且标注上了样本的方差：

1
2
3
4
5
6
7
8
9
10
11
12
library(rafalib)
mypar(3,3)
sds=seq(0,2,len=100)
for(d in c(1,5,10)){
  for(s0 in c(0.1, 0.2, 0.3)){
    tmp=hist(biosds,main=paste("s_0 =",s0,"d =",d),xlab="sd",ylab="density",freq = FALSE,nc=100,xlim=c(0,1))
    dd=df(sds^2/s0^2,11,d)
    ##multiply by normalizing constant to assure same range on plot
    k=sum(tmp$density)/sum(dd)
    lines(sds,dd*k,type="l",col=2,lwd=2)
  }
}

从上面的图形我们要找到哪个图形的$s_{0}$与$d$最能拟合我们的数据？这是一个比较麻烦的过程，因为MLE并不会计算这个特殊的分布（具体可以参考Wmyth(2004)）。Bioconductor上的limma包提供的估计这些参数的函数，如下所示：

1
2
3
4
library(limma)
estimates=fitFDist(biosds^2,11)
theoretical<- sqrt(qf((seq(0,999)+0.5)/1000, 11, estimates$df2)*estimates$scale)
observed <- biosds

拟合模型确实提供了一个合理的估计，从qq图和直方图就能看出来：

1
2
3
4
5
6
7
mypar(1,2)
qqplot(theoretical,observed)
abline(0,1)
tmp=hist(biosds,main=paste("s_0 =", signif(estimates[[1]],2), "d =", signif(estimates[[2]],2)), xlab="sd", ylab="density", freq=FALSE, nc=100, xlim=c(0,1), ylim=c(0,9))
dd=df(sds^2/estimates$scale,11,estimates$df2)
k=sum(tmp$density)/sum(dd) ##a normalizing constant to assure same area in plot
lines(sds, dd*k, type="l", col=2, lwd=2)

练习

P297